library(specr)
library(fixest)
library(readr)
library(broom.mixed)
library(texreg)
library(dotwhisker)
library(dplyr)
library(fixest)

# Two-Way logistic Fixed Effect Models that use all data
# but due to logit characteristics drop observations where there is no variation
# across time in the outcomes

set.seed(42)

source("funs_and_constants.R")

# add in data
epr_yearly_group_clusters_h1 <- read_csv("data/epro_data_h1.csv")
epro_year_group_clusters <- read_csv("data/epro_data_h2.csv")
epro_year_group_clusters_geo <- read_csv("data/epro_geo_data_h3.csv")

full_model_h1 <- feglm(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +  n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                       data = epr_yearly_group_clusters_h1, family = "binomial")

full_model_h1 <- feglm(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +  n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                       data = epr_yearly_group_clusters_h1, family = "binomial")

full_model_h1_n_resources <- feglm(disagreement ~ none_one_or_both_nats_agri_char + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                                   +  n_tek_groups | year+ countries_gwid, cluster = ~countries_gwid, data = epr_yearly_group_clusters_h1, family = "binomial")


model_list_h1 <- list(full_model_h1, full_model_h1_n_resources)

texreg::screenreg(model_list_h1, stars = stars, custom.coef.map = coef_map)

texreg(model_list_h1,
       stars = stars,
       custom.coef.map = coef_map,
       table = FALSE,
       custom.model.names = c("Disagreement (1/0)", "Disagreement (1/0)"),
       include.proj.stats = FALSE,
       file = "tables/h1_logit_panel.tex")


##########
##### H2
##########

full_model_h2 <- feglm(disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +  n_tek_groups | year +countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters, family = "binomial")

# n rels as IV
full_model_h2_n_rels <- feglm(disagreement ~ n_ed_religions + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                              +  n_tek_groups | year + countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters, family = "binomial")

full_model_h2_herf <- feglm(disagreement ~ hhi_rel + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                            +  n_tek_groups | year + countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters, family = "binomial")


screenreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf),
          stars = stars)

texreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:3),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h2_logit_panel.tex")

######
### h3
#######

full_model_h3 <- feglm(disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +  n_tek_groups | year + countries_gwid,
                       cluster = ~countries_gwid,
                       data = epro_year_group_clusters_geo,
                       family = "binomial"
)

full_model_h3_bin <- feglm(disagreement ~ multiple_polygons_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                           +  n_tek_groups | year + countries_gwid,
                           cluster = ~countries_gwid,
                           data = epro_year_group_clusters_geo,
                           family = "binomial"
)

full_model_h3_hhi <- feglm(disagreement ~ spatial_hhi + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                           +  n_tek_groups | year + countries_gwid,
                           cluster = ~countries_gwid,
                           data = epro_year_group_clusters_geo,
                           family = "binomial"
)



texreg(list(full_model_h3, full_model_h3_bin, full_model_h3_hhi),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:3),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h3_logit_panel.tex")





